Użyte biblioteki

library(readxl)
library(httr)
library(arules)
library(tidyr)
library(dplyr)
library(lubridate)
library(stringr)
library(ggplot2)
library(purrr)
library(broom)
library(tibble)
library(plotly)

Przetwarzanie zbioru danych

Odczyt danych

Odczyt danych przeprowdzony jest bezośrednio ze strony, aby zapewenić stałą wartość początkową dla każdego uruchomienia. Dodatkowo uzupełniane są brakujące dane identyfikatory pacjentów PATIENT_ID. Wartości brakujące wynikają ze sposobu przetwarzania scalonych komórek pliku źródłowego.

# https://stackoverflow.com/a/41368947
url <- "http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx"
GET(url, write_disk(tf <- tempfile(fileext = ".xlsx")))
df <- read_excel(tf) %>% fill(PATIENT_ID)

Atrybuty zbioru danych

Surowy zbiór danych zawiera 81 atrybutów i 6120 wierszy.

Atrybuty zbioru danych można potwierdzić na dwie grupy:

  • Atrybuty ogólne - zawierają etykiety, podstawowe dane pacjenta(płeć, wiek, identyfikator, datę przyjęcia i wypisu, rezultat) oraz datę wykonania badania
  • Wyniki badań - w kolumnach są umieszczone wartości poszczególnych badanych właściwości krwi

Brakujące dane

Poza danymi niedostępnymi (NA) w zbiorze danych wsytępują dane brakujące, oznakowane wartością -1. Poniższa tabela pokazuje ilość takich wartości w każdej kolumnie, o ile te wartości występują.

df %>% pivot_longer(-c(RE_DATE,`Admission time`, `Discharge time`)) %>% filter(value==-1) %>% count(name)
## # A tibble: 2 x 2
##   name                                 n
##   <chr>                            <int>
## 1 2019-nCoV nucleic acid detection   501
## 2 Platelet count                       4

Dla kolumny Platelet count warkości brakujących jest mało, liczby -1 zostaną zamienione na NA.

df <- df %>% mutate(`Platelet count` = if(is.na(`Platelet count`) || `Platelet count` != -1) `Platelet count` else NA)

Dla kolumny 2019-nCoV nucleic acid detection sprawdzimy całkowitą ilość wartości różnych od NA.

df %>% 
  pivot_longer(-c(RE_DATE,`Admission time`, `Discharge time`)) %>% 
  filter(name=="2019-nCoV nucleic acid detection") %>% 
  filter(!is.na(value)) %>% count(name)
## # A tibble: 1 x 2
##   name                                 n
##   <chr>                            <int>
## 1 2019-nCoV nucleic acid detection   501

Wszystkie wartości w kolumnie 2019-nCoV nucleic acid detection różne od NA wynoszą -1, w związku z tym ta kolumna jest właściwie nieprzydatna. Mimo to na ten moment postanowiono nie zmieniać wartości w tej kolumnie, ewentualne zmiany można wykonać w dalszych krokach.

Atrybuty ogólne

Pierwszych siedem kolumn to atrybuty ogólne, reprezentują następująco:

  • PATEINT_ID - sztuczny identyfikator pacjenta
  • RE_DATE - najpewniej jest to data wprowadzenia wyniku badań próbki krwi do systemu
  • age - wiek pacjenta
  • gender - płeć pacjenta, przyjmuje wartość 1 albo 2, ale z poisu zbioru danych nie wynika wprost, które oznaczenie jest przypisane do której płci
  • Admission time - czas przyjęcia pacjenta do szpitala
  • Discharge time - czas zakończenia leczenia (w skutek wypisu pacjenta lub jego śmierci)
  • outcome - rezultat, przyjmuje wartość 0 jeżeli pacjent został wyleczony, lub wartość 1 jeżeli zmarł.

Nazwy wyżej wymienionych atrybutów zostaną przeformatowane, będą pisane wielkimi literami, a ewentualne spacje będą zamieniane na znak podkreślinka _, w ten sposób atrbuty ogólne będą łatwo odróżnialne od pozostałych atrybutów.

df <- df %>% rename_with(~ toupper(gsub(" ", "_", .x, fixed = TRUE)), 1:7)

Pola GENDER i OUTCOME zostaną dodatkowo zmienione na typ factor - z typu numerycznego zamienione zostaną na zmienne nominalne. Do tego trzeba przeprowadzić identyfikację typu płci.

Wartości płci nie są jawnie opisane, w artykule źródłowym została natomiast podana proporcja pacjentów danej płci, co powinno umożliwić okreslenie właściwego przypisania.

knitr::kable(
df %>% group_by(PATIENT_ID) %>% summarise_all(last) %>% count(GENDER) %>% rename("count" = "n"),
align="lc"
)
GENDER count
1 224
2 151

Wniosek:

  • Mężczyźni - wartość 1
  • Kobiety - wartość 2
df <- df %>% mutate(across("GENDER", ~factor(., levels=c(1,2), labels=c("MALE", "FEMALE")))) %>%
  mutate(across("OUTCOME", ~factor(., levels=c(0,1), labels=c("CURED","DECEASED"))))

Poniżej można zapoznać się z analizą wartości dla atrybutów ogólnych. Wartość atrybutu RE_DATE ma unikalną wartość dla kazdego wiersza.

knitr::kable(summary(select(df, RE_DATE)))
RE_DATE
Min. :2020-01-10 19:45:00
1st Qu.:2020-02-04 13:44:00
Median :2020-02-09 12:42:30
Mean :2020-02-08 07:00:02
3rd Qu.:2020-02-13 10:34:00
Max. :2020-02-18 17:49:00
NA’s :14

Pozostałe z atrybtów ogólnych, są wspólne dla każdego pacjenta, zastosowano transformację do zniwelowania skutków duplikacji.

options(knitr.kable.NA = '')
knitr::kable(
  summary(
    df %>% 
      select(c(1,3:7)) %>% 
      group_by(PATIENT_ID) %>% 
      summarise_all(last)
    )
  )
PATIENT_ID AGE GENDER ADMISSION_TIME DISCHARGE_TIME OUTCOME
Min. : 1.0 Min. :18.00 MALE :224 Min. :2020-01-10 15:52:20 Min. :2020-01-23 09:09:23 CURED :201
1st Qu.: 94.5 1st Qu.:46.00 FEMALE:151 1st Qu.:2020-02-01 19:27:40 1st Qu.:2020-02-11 13:39:21 DECEASED:174
Median :188.0 Median :62.00 Median :2020-02-04 22:30:34 Median :2020-02-16 17:40:07
Mean :188.0 Mean :58.83 Mean :2020-02-04 20:13:51 Mean :2020-02-15 16:42:59
3rd Qu.:281.5 3rd Qu.:70.00 3rd Qu.:2020-02-10 04:11:10 3rd Qu.:2020-02-19 11:47:14
Max. :375.0 Max. :95.00 Max. :2020-02-17 21:30:07 Max. :2020-03-04 16:21:51

Liczba unikatowych identyfikatorów pacjentów wynosi 375.

Czas hosptializacji

Ta sekcja poświęcona jest analizie długości hospitalizacji poszczególnych grup pacjentów.

HOSPITALIZATION_TIME quantile
1 days min
6 days 1st
11 days 2nd
17 days 3rd
36 days max

OUTCOME GENDER AGE HOSPITALIZATION_TIME
CURED MALE 54.5 15.0 days
CURED FEMALE 45.0 15.0 days
DECEASED MALE 68.0 7.0 days
DECEASED FEMALE 72.0 6.5 days

Wyniki badań

Pozostałe 74 kolumn zawiera wartości atrybutów, które reprezentują rodzaje miar zebranych podczas każdego badania.

Nazwy kolumn zostały poddane następującym transformacjom:

  • Zamieniono długie nazwy kolumn na krótkie, skrótowe odpowiedniki, na podstawie nazw medycznych (częściowy code book dostępny jest w repozytorium)
  • oznaczenia (%) zamieniono na suffix _percent
  • oznaczenia (#) usunięto
  • nazwy zaczynajace się wielką literą i poza tym skłądające się z małych liter zamieniono na odpowiedniki skłądające się wyłącznie z małych znaków, jednak skrótowce zostały zachwowane
  • znak spacji i myślnika - zamieniono na podkreślenie dolne _ ze względu na łatwość odwołania się do kolumn

Poniższy blok kodu szczegółowo ukazuje zatosowane transfomrmacje, jednocześnie prezentowana są ostatenczne nazwy kolumn wyników badań.

df <- df %>%
  rename("NT-proBNP" = "Amino-terminal brain natriuretic peptide precursor(NT-proBNP)") %>%
  rename("2019-nCov-detection" = "2019-nCoV nucleic acid detection") %>%
  rename("aPPT" = "Activation of partial thromboplastin time") %>%
  rename("ALP" = "Alkaline phosphatase") %>%
  rename("AAT" = "aspartate aminotransferase") %>%
  rename("FDPs" = "Fibrin degradation products") %>%
  rename("ALT" = "glutamic-pyruvic transaminase") %>%
  rename("hs-CRP" = "High sensitivity C-reactive protein") %>%
  rename("hs-cTnI" = "Hypersensitive cardiac troponinI") %>%
  rename("HCV-abs-quant" = "HCV antibody quantification") %>%
  rename("HIV-abs-quant" = "HIV antibody quantification") %>%
  rename("INR" = "International standard ratio") %>%
  rename("LDH" = "Lactate dehydrogenase") %>%
  rename("MCH" = "mean corpuscular hemoglobin") %>%
  rename("MCHC" = "mean corpuscular hemoglobin concentration") %>%
  rename("MCV" = "mean corpuscular volume") %>%
  rename("MPV" = "Mean platelet volume") %>%
  rename("P-LCR" = "platelet large cell ratio") %>%
  rename("PDW" = "PLT distribution width") %>%
  rename("q-t-pallidum-abs" = "Quantification of Treponema pallidum antibodies") %>%
  rename("RCDW-SD" = "RBC distribution width SD") %>%
  rename("RBC_count" = "Red blood cell count") %>%
  rename("RCDW" = "Red blood cell distribution width") %>%
  rename("TNF-alfa" = "Tumor necrosis factorα") %>%
  rename("WBC_count" = "White blood cell count") %>%
  rename("gamma-GT" = "γ-glutamyl transpeptidase") %>%
  rename_with(~str_c(str_replace(.x, "\\(%\\)", ""), "_percent"), contains("(%)")) %>%
  rename_with(~str_replace(.x, "\\(#\\)", "")) %>%
  rename_with(tolower, matches("^[[:upper:]]{1}[[:lower:]]+", ignore.case = FALSE)) %>%
  rename_with(~str_replace_all(.x, "[ |-]", "_"))
str_sort(colnames(df[,-(1:7)]))
##  [1] "2019_nCov_detection"    "AAT"                    "albumin"               
##  [4] "ALP"                    "ALT"                    "antithrombin"          
##  [7] "aPPT"                   "basophil_count"         "basophil_percent"      
## [10] "calcium"                "corrected_calcium"      "creatinine"            
## [13] "D_D_dimer"              "direct_bilirubin"       "eGFR"                  
## [16] "eosinophil_count"       "eosinophils_percent"    "ESR"                   
## [19] "FDPs"                   "ferritin"               "fibrinogen"            
## [22] "gamma_GT"               "globulin"               "glucose"               
## [25] "HBsAg"                  "HCO3_"                  "HCV_abs_quant"         
## [28] "hematocrit"             "hemoglobin"             "HIV_abs_quant"         
## [31] "hs_CRP"                 "hs_cTnI"                "indirect_bilirubin"    
## [34] "INR"                    "interleukin_10"         "interleukin_1β"        
## [37] "interleukin_2_receptor" "interleukin_6"          "interleukin_8"         
## [40] "LDH"                    "lymphocyte_count"       "lymphocyte_percent"    
## [43] "MCH"                    "MCHC"                   "MCV"                   
## [46] "monocytes_count"        "monocytes_percent"      "MPV"                   
## [49] "neutrophils_count"      "neutrophils_percent"    "NT_proBNP"             
## [52] "P_LCR"                  "PDW"                    "PH_value"              
## [55] "platelet_count"         "procalcitonin"          "prothrombin_activity"  
## [58] "prothrombin_time"       "q_t_pallidum_abs"       "RBC_count"             
## [61] "RCDW"                   "RCDW_SD"                "serum_chloride"        
## [64] "serum_potassium"        "serum_sodium"           "thrombin_time"         
## [67] "thrombocytocrit"        "TNF_alfa"               "total_bilirubin"       
## [70] "total_cholesterol"      "total_protein"          "urea"                  
## [73] "uric_acid"              "WBC_count"

Czyszczenie zbioru danych

Niniejsza sekcja jest poświęcona wstępnemu przetwarzniu zbiorów danych, ze szczególnym uwzględnieniem usuwania i scalania wybranych rekordów.

Puste rekordy w zbiorze

W powyższej głównej zauważono wystąpienia wartości brakujących dla kolumny RE_DATE. Szybkie sprawdzenie wykazuje, że dla wierszy, gdzie nie ma atrybut RE_DATE przyjmuje wartość ‘NA’, wartości wszystkich atrybutów pochodzących z badań krwi, również przyjmują wartość ‘NA’.

all(sapply(df[is.na(df$RE_DATE), -c(1,3:7)], is.na))
## [1] TRUE

Fakt, że wiersze te są właściwie puste, jest podstawą do ich usunięcia ze zbioru.

df <- df[!is.na(df$RE_DATE),]

Po tej zmianie w zbiorze pozostało 6106 wierszy i 361 unikalnych identyfikatorów pacjentów.

Kategorie badań

Wstępnę analiza danych sugeruje wskazuje, że wiele wartości atrybutów jest oznaczona wartością ‘NA’. Poniższy wykresy wskazuje, ilość wystąpień rekordów badań z danę liczbą wartości (różnych od NA).

Z powyższego histogramu można wyciągnąć kilka wniosków

  • Znacząca większość rekordów ma zaledwie jeden wynik badania (jedną wartość).
  • Interesujący jest fakt, że druga i trzecia pod względem liczości grupa rekordów ma odpowiednio po 24 i 23 elementy, taki fakt może sugerować, że elementy mogą występować w grupach. Możliwa jest dalsza analiza współwystępowania ze sobą atrybutów, np. przy pomocy analizy zbiorów częstych.
  • Ze względu na to, że właściwie nie ma rekordów, które opisują wszystkie 74 atrybuty, konieczne będzie scalanie rekordów w dalszych krokach.

Współwystępowanie kategorii

W poprzedniej sekcji stwierdzono, że można przeanalizować współwystępowanie ze sobą poszczególnych atrybutów w rekordach. W związku z tym tabela wyników badań zostanie przekształcona do zbioru binarnych transakcji.

transactions <- sapply(df[,-(1:7)], function(x) !is.na(x))

Sprawdzona zostanie ilość unikalnych transakcji.

length(unique(apply(transactions, 1, function(x) which(x))))
## [1] 176

W zbiorze 6106 transakcji jest zaledwie 176 różnych typów.

Aby sprawdzić, czy atrybuty wsytępują w grupach, wyszukiwane są maksymalne domknięte zbiory częste, dla mimalnej wartości support równej 5%.

itemsets <- apriori(transactions, parameter = list(target="closed frequent itemsets", support=.05, minlen=1, maxlen=30, maxtime=60))
## Apriori
## 
## Parameter specification:
##  confidence minval smax arem  aval originalSupport maxtime support minlen
##          NA    0.1    1 none FALSE            TRUE      60    0.05      1
##  maxlen                   target  ext
##      30 closed frequent itemsets TRUE
## 
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
## 
## Absolute minimum support count: 305 
## 
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[74 item(s), 6106 transaction(s)] done [0.00s].
## sorting and recoding items ... [63 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 done [26.87s].
## filtering closed item sets ... done [15.64s].
## sorting transactions ... done [0.00s].
## writing ... [130 set(s)] done [0.34s].
## creating S4 object  ... done [1.01s].
mc_itemsets <- itemsets[is.maximal(itemsets)]
inspect(sort(mc_itemsets, by = 'support'))
##      items                     support transIdenticalToItemsets count
## [1]  {hemoglobin,                                                    
##       eosinophils_percent,                                           
##       basophil_percent,                                              
##       platelet_count,                                                
##       monocytes_percent,                                             
##       RCDW,                                                          
##       neutrophils_percent,                                           
##       MCV,                                                           
##       hematocrit,                                                    
##       WBC_count,                                                     
##       MCHC,                                                          
##       lymphocyte_count,                                              
##       RBC_count,                                                     
##       eosinophil_count,                                              
##       neutrophils_count,                                             
##       MPV,                                                           
##       RCDW_SD,                                                       
##       lymphocyte_percent,                                            
##       P_LCR,                                                         
##       monocytes_count,                                               
##       PDW,                                                           
##       basophil_count,                                                
##       MCH,                                                           
##       thrombocytocrit}      0.13609564                0.1354406   831
## [2]  {glucose}              0.12692434                0.0000000   775
## [3]  {serum_chloride,                                                
##       ALP,                                                           
##       albumin,                                                       
##       total_bilirubin,                                               
##       indirect_bilirubin,                                            
##       total_protein,                                                 
##       urea,                                                          
##       corrected_calcium,                                             
##       serum_potassium,                                               
##       direct_bilirubin,                                              
##       total_cholesterol,                                             
##       AAT,                                                           
##       uric_acid,                                                     
##       HCO3_,                                                         
##       calcium,                                                       
##       LDH,                                                           
##       globulin,                                                      
##       gamma_GT,                                                      
##       hs_CRP,                                                        
##       serum_sodium,                                                  
##       ALT,                                                           
##       eGFR,                                                          
##       creatinine}           0.11021946                0.1094006   673
## [4]  {hs_cTnI}              0.08303308                0.0000000   507
## [5]  {2019_nCov_detection}  0.08205044                0.0000000   501
## [6]  {NT_proBNP}            0.07779234                0.0000000   475
## [7]  {procalcitonin}        0.07517196                0.0000000   459
## [8]  {PH_value}             0.06288896                0.0000000   384
## [9]  {ESR}                  0.06272519                0.0000000   383
## [10] {prothrombin_time,                                              
##       antithrombin,                                                  
##       prothrombin_activity,                                          
##       fibrinogen,                                                    
##       thrombin_time,                                                 
##       D_D_dimer,                                                     
##       FDPs,                                                          
##       INR,                                                           
##       aPPT}                 0.05371765                0.0000000   328

Otrzymaliśmy 10 maksymalnych domkniętych zbiorów częstych, z czego 7 to zbiory jednolementowe. Pozostałe 3 zbiory (1, 3 i 10) zawierają wiele elementów. Na podstawie elementów składowych, zbiory wieolemelemnetowe zostały zakwalifikowane następująco:

  • Zbiór 1 - zawiera elementy typowe dla badania morfologii krwi (ang. Compete blood count (CBC))
  • Zbiór 3 - zawiera elementy typowe dla badania biochemii krwi (ang. Blood biochemistry)
  • Zbiór 10 - zawiera elementy typowe dla badania krzepliwości krwi (ang. Coagulation)

Na podstawie uzyskanego wyniku można wnioskować, że ilość elementów w rekordzie może wynikać z procesu technologicznego analizowania próbek krwi. Sugeruje to, że możliwe będzie scalanie rekordów, np na podstawie dnia wykonania badania (jeśli jednemu pacjentowi przeprowadzono wiele badań w ciagu dnia).

Badania w czasie

TODO

class_labels <- as(items(sort(mc_itemsets, by="support")), "list")
class_labels <- sapply(class_labels, function(x) x[1])
class_labels[1] <- "Complete blood count"
class_labels[3] <- "Blood biochemistry"
class_labels[10] <- "Coagulation"
classes <- as(items(sort(mc_itemsets, by="support")), "matrix")
labeled <- df %>% rowwise() %>%
  mutate(CLASS = paste(
    which(
      apply(
        (matrix(
          !is.na(c_across(-(1:7))), 
          nrow=nrow(classes), 
          ncol=ncol(classes), 
          byrow=TRUE, 
          dimnames=dimnames(classes)
          ) & classes) == classes, 
        1, 
        all)
      ), collapse=""), 
    .after=OUTCOME ) %>%
  mutate(CLASS = if(CLASS != "" && as.integer(CLASS)>10) "MULTIPLE" else CLASS) %>%
  mutate(CLASS = if(CLASS != "" && CLASS != "MULTIPLE") class_labels[as.integer(CLASS)] else CLASS) %>%
  mutate(CLASS = str_replace(CLASS, "^$", "NONE")) %>%
  mutate(CLASS = str_trim(CLASS)) %>%
  mutate(CLASS = factor(CLASS)) %>%
  rowwise() %>%
  mutate(SIZE = sum(!is.na(c_across(-c(1:8)))), .after=CLASS) %>%
  select(1:9) %>%
  ungroup()

labeled %>% group_by(CLASS) %>% summarize(CNT = n(), .groups="drop") %>% arrange(desc(CNT))
## # A tibble: 12 x 2
##    CLASS                  CNT
##    <fct>                <int>
##  1 NONE                  1433
##  2 Complete blood count   816
##  3 glucose                584
##  4 Blood biochemistry     507
##  5 2019_nCov_detection    500
##  6 MULTIPLE               485
##  7 hs_cTnI                386
##  8 PH_value               375
##  9 ESR                    368
## 10 Coagulation            300
## 11 NT_proBNP              190
## 12 procalcitonin          162
group.plot <- ggplot(
  labeled, 
  aes(x=RE_DATE, y=CLASS, color=CLASS)
  ) +
  geom_jitter(size=0.2) +
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b %Y") +
  theme(legend.position="none")

ggplotly(group.plot)

Badania poszczególnych pacjentów

TODO Analiza ilości badań przeprowadzonych jednemu pacjentowi.

patient_last_results <- df %>% 
  group_by(PATIENT_ID) %>% 
  fill(everything()) %>% 
  summarise_all(last) %>%
  rowwise() %>% 
  mutate(number_of_tests=sum(!is.na(c_across(-(1:7)))), .after=OUTCOME)
  
ggplot(
  patient_last_results, 
  aes(number_of_tests)
  ) + geom_histogram(binwidth = 1)

Zauważamy, że część pacjentów ma niewiele przeprowadzonych rodzajów badań.

TODO analityczne wybranie tej wartości
* procent pełnych kolumn * procent pozostawionych pacjentów * wartość progu

veryfication_df <- function(n_tests, cols_cnt, th_seq) {
  cls_vec <- c()
  pts_vec <- c()
  for(th in th_seq) {
    tmp <- n_tests %>% filter(number_of_tests < cols_cnt*th)
    dt <- n_tests %>% filter(!(PATIENT_ID %in% tmp$PATIENT_ID))
    fl_cols <- mean(colSums(sapply(dt[,-(1:8)], is.na)) == 0)
    pnts <- length(unique(dt$PATIENT_ID)) / length(unique(n_tests$PATIENT_ID))
    cls_vec <- append(cls_vec, fl_cols)
    pts_vec <- append(pts_vec, pnts)
  }
    data.frame(threshold=th_seq, full_cols=cls_vec, patients_left=pts_vec)
}

to_plot <- veryfication_df(patient_last_results, ncol(df)-8, seq(0.5, 1.0, by=.025))

ggplot(to_plot, aes(x=threshold)) +
  geom_line(aes(y=full_cols, color="Full columns")) + 
  geom_line(aes(y=patients_left, color="Patients left")) + 
  labs(color="Legend") +
  geom_vline(xintercept=0.675, linetype = "longdash")

TODO

W zbiorze danych pozostawimy tylko tych, którzy mają ponad 0.675 wykonanych rodzajów badań.

threshold <- 0.675
patients_to_be_removed <- patient_last_results %>% 
  filter(number_of_tests < (ncol(df)-7) * threshold) %>% 
  select(PATIENT_ID)

df <- df %>% filter(PATIENT_ID %!in% patients_to_be_removed$PATIENT_ID)

Po usunięciu rekordów pacjentów, u których przeprowadzono mało badań, pozostało 354 unikalnych pacjentów.

Brakujące wartości w kategoriach

knitr::kable(
df %>% 
  group_by(PATIENT_ID) %>% 
  fill(everything()) %>% 
  summarise_all(last) %>% 
  ungroup() %>%
  select(-(1:7)) %>% 
  pivot_longer(everything()) %>%
  mutate(value=is.na(value)) %>%
  group_by(name) %>%
  summarize(value=sum(value), .groups="drop") %>%
  arrange(value)
)
name value
AAT 0
albumin 0
ALP 0
ALT 0
basophil_count 0
basophil_percent 0
creatinine 0
direct_bilirubin 0
eGFR 0
eosinophil_count 0
eosinophils_percent 0
gamma_GT 0
globulin 0
HCO3_ 0
hematocrit 0
hemoglobin 0
LDH 0
lymphocyte_count 0
lymphocyte_percent 0
MCH 0
MCHC 0
MCV 0
monocytes_count 0
monocytes_percent 0
neutrophils_count 0
neutrophils_percent 0
platelet_count 0
RBC_count 0
total_bilirubin 0
total_cholesterol 0
total_protein 0
urea 0
uric_acid 0
WBC_count 0
indirect_bilirubin 1
calcium 2
serum_chloride 2
serum_potassium 2
serum_sodium 2
corrected_calcium 3
hs_CRP 3
glucose 5
INR 5
prothrombin_activity 5
prothrombin_time 5
RCDW 6
RCDW_SD 6
MPV 10
P_LCR 10
PDW 10
thrombocytocrit 10
D_D_dimer 14
procalcitonin 44
aPPT 59
fibrinogen 59
thrombin_time 59
ESR 69
hs_cTnI 69
HBsAg 82
HCV_abs_quant 82
q_t_pallidum_abs 82
HIV_abs_quant 83
NT_proBNP 90
PH_value 126
interleukin_6 138
2019_nCov_detection 140
interleukin_1β 140
interleukin_2_receptor 140
interleukin_8 140
TNF_alfa 140
interleukin_10 141
ferritin 143
antithrombin 153
FDPs 153

Dodatkowe cechy

  • TODO days to outcome
df <- df %>%  
  mutate(DAYS_TO_OUTCOME = as.numeric(difftime(
    floor_date(DISCHARGE_TIME, "1 day"), 
    floor_date(RE_DATE, "1 day"), unit="days"
    ), unit="days"), .before=OUTCOME)
df %>% summarise(DAYS_TO_OUTCOME = quantile(DAYS_TO_OUTCOME, c(0, 0.25, 0.5, 0.75,1)), quantile = c("min", "1st", "2nd", "3rd", "max"))
## # A tibble: 5 x 2
##   DAYS_TO_OUTCOME quantile
##             <dbl> <chr>   
## 1              -2 min     
## 2               3 1st     
## 3               7 2nd     
## 4              13 3rd     
## 5              35 max
ggplot(df, aes(DAYS_TO_OUTCOME)) + geom_histogram(binwidth = 1)

Agregacja wierszy według daty badania

df <- df %>% 
  mutate(RE_DATE = floor_date(RE_DATE, '1 day')) %>% 
  group_by(PATIENT_ID, RE_DATE) %>% 
  fill(everything()) %>% 
  summarise_all(last) %>% 
  ungroup() %>%
  rowwise() %>%
  mutate(number_of_tests=sum(!is.na(c_across(-(1:8)))), .after=OUTCOME)
ggplot(df, aes(number_of_tests)) + geom_histogram(binwidth = 1)

df <- df %>% select(-number_of_tests)

Zbiór danych

ggplot(
  df %>% pivot_longer(hs_cTnI:creatinine) %>% filter(!is.na(value)),
  aes(x=OUTCOME, y=value)
) + geom_boxplot() +
scale_x_discrete(labels = abbreviate) + 
facet_wrap(name ~ ., scales="free",ncol=5)

Korelacja

TODO

cor_mat <- cor(
  x = df %>% 
    mutate(isMALE = as.numeric(GENDER=="MALE"), 
           isCURED = as.numeric(OUTCOME=="CURED"),
           has_nCOV_test = as.numeric(!is.na(`2019_nCov_detection`))) %>% 
    select(-c(1:8), -"2019_nCov_detection"), 
  use="pairwise.complete.obs"
  )
cor_df = data.frame(round(cor_mat,2)) %>% 
 rownames_to_column() %>%
 pivot_longer(-rowname, names_to="colname")

cor_plot <- ggplot(cor_df, aes(colname, rowname, fill=value)) + 
  geom_tile() +
  scale_fill_gradient2() + 
  theme(axis.text.x=element_text(angle = 90, hjust = 0))

ggplotly(cor_plot)
knitr::kable(
cor_df %>% filter(colname>rowname) %>% arrange(desc(abs(value))) %>% head(30)
)
rowname colname value
MPV P_LCR 0.99
INR prothrombin_time 0.99
platelet_count thrombocytocrit 0.98
HBsAg interleukin_6 0.98
direct_bilirubin total_bilirubin 0.98
lymphocyte_percent neutrophils_percent -0.98
procalcitonin q_t_pallidum_abs 0.95
D_D_dimer FDPs 0.95
P_LCR PDW 0.95
MPV PDW 0.94
ALT interleukin_1β 0.93
lymphocyte_count monocytes_count 0.88
serum_chloride serum_sodium 0.86
eosinophil_count eosinophils_percent 0.86
AAT interleukin_1β 0.86
hematocrit hemoglobin 0.84
MCH MCV 0.82
indirect_bilirubin total_bilirubin 0.80
interleukin_6 interleukin_8 0.80
RCDW RCDW_SD 0.79
aPPT prothrombin_time 0.79
monocytes_percent neutrophils_percent -0.77
eGFR urea -0.77
creatinine urea 0.75
albumin calcium 0.74
isCURED neutrophils_percent -0.73
isCURED lymphocyte_percent 0.72
neutrophils_count neutrophils_percent 0.71
albumin total_protein 0.70
lymphocyte_percent neutrophils_count -0.70
knitr::kable(
tidy(cor.test(df$LDH, as.numeric(df$OUTCOME=="CURED")))
)
estimate statistic p.value parameter conf.low conf.high method alternative
-0.6270327 -24.33491 0 914 -0.6648056 -0.5860615 Pearson’s product-moment correlation two.sided
correlations <- df %>% mutate(isCURED=as.numeric(OUTCOME=="CURED")) %>% select(-c(1:2,4:8)) %>% 
  pivot_longer(!isCURED, names_to="attribute", values_to="x_val") %>%
  filter(!is.na(x_val)) %>%
  nest(data = c(x_val, isCURED)) %>%
  mutate(cor_test = map(data, ~tidy(cor.test(.x$x_val, .x$isCURED)))) %>% 
  unnest(cor_test) %>%
  select(attribute, estimate, p.value, conf.low, conf.high) %>% 
  arrange(desc(abs(estimate))) %>% head(15)
## Warning: Problem with `mutate()` input `cor_test`.
## ℹ the standard deviation is zero
## ℹ Input `cor_test` is `map(data, ~tidy(cor.test(.x$x_val, .x$isCURED)))`.
## Warning in cor(x, y): the standard deviation is zero
knitr::kable(
correlations
)
attribute estimate p.value conf.low conf.high
neutrophils_percent -0.7301461 0 -0.7586392 -0.6988654
lymphocyte_percent 0.7234883 0 0.6915912 0.7525691
albumin 0.6880859 0 0.6524294 0.7207028
prothrombin_activity 0.6547754 0 0.6084767 0.6966318
hs_CRP -0.6509896 0 -0.6910468 -0.6069460
D_D_dimer -0.6376192 0 -0.6819778 -0.5885867
LDH -0.6270327 0 -0.6648056 -0.5860615
neutrophils_count -0.6083675 0 -0.6470960 -0.5665074
FDPs -0.6042409 0 -0.6689583 -0.5304310
AGE -0.5761382 0 -0.6071595 -0.5433633
calcium 0.5703332 0 0.5257539 0.6117881
platelet_count 0.5419441 0 0.4952125 0.5855487
monocytes_percent 0.5415337 0 0.4947996 0.5851444
eGFR 0.4909936 0 0.4402769 0.5385870
urea -0.4841428 0 -0.5321758 -0.4330031

Analiza wypisów ze szpitala

tseries <- df %>% 
  select(1:8) %>% 
  group_by(PATIENT_ID) %>% 
  summarize(across(everything(), last), .groups="drop") %>% 
  select(-RE_DATE, -AGE, -GENDER, -DAYS_TO_OUTCOME, -PATIENT_ID) %>%
  pivot_longer(ADMISSION_TIME:DISCHARGE_TIME, names_to="TYPE", values_to="DATE") %>%
  mutate(DATE = floor_date(DATE, "days")) %>%
  mutate(TYPE = str_replace(TYPE, "_TIME", "")) %>% 
  count(TYPE, DATE, OUTCOME) %>%
  rename("COUNT" = "n")

tplot <- ggplot(tseries, aes(x=DATE, y=COUNT, fill=OUTCOME)) + 
  geom_col() + 
  scale_x_datetime(date_breaks = "1 week", date_labels = "%d %b %Y") +
  facet_grid(TYPE ~ .)

ggplotly(tplot)

Techniczne

Wszytkie użyte pakiety

unique(c(.packages(), loadedNamespaces()))
##  [1] "plotly"      "tibble"      "broom"       "purrr"       "ggplot2"    
##  [6] "stringr"     "lubridate"   "dplyr"       "tidyr"       "arules"     
## [11] "Matrix"      "httr"        "readxl"      "stats"       "graphics"   
## [16] "grDevices"   "utils"       "datasets"    "methods"     "base"       
## [21] "Rcpp"        "highr"       "cellranger"  "compiler"    "pillar"     
## [26] "tools"       "digest"      "viridisLite" "jsonlite"    "evaluate"   
## [31] "lifecycle"   "gtable"      "lattice"     "pkgconfig"   "rlang"      
## [36] "cli"         "crosstalk"   "curl"        "yaml"        "xfun"       
## [41] "withr"       "knitr"       "htmlwidgets" "generics"    "vctrs"      
## [46] "grid"        "tidyselect"  "data.table"  "glue"        "R6"         
## [51] "fansi"       "rmarkdown"   "farver"      "magrittr"    "backports"  
## [56] "scales"      "htmltools"   "ellipsis"    "assertthat"  "colorspace" 
## [61] "labeling"    "utf8"        "stringi"     "lazyeval"    "munsell"    
## [66] "crayon"